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Abstract 

We review recent dissipative particle dynamics (DPD) simulations of electrolyte flow in nanochannels. A method is 
presented by which the slip length at the channel boundaries can be tuned systematically from negative to infinity 
by introducing suitably adjusted wall-fluid friction forces. Using this method, we study electroosmotic flow (EOF) in 
nanochannels for varying surface slip conditions and fluids of different ionic strength. Analytic expressions for the 
flow profiles are derived from the Stokes equation, which are in good agreement with the numerical results. Finally, 
we investigate the influence of EOF on the effective mobility of polyelectrolytes in nanochannels. The relevant 
quantity characterizing the effect of slippage is found to be the dimensionless quantity k6b, where l//c is an effective 
electrostatic screening length at the channel boundaries. 
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1. Introduction 

Microfluidic devices like bio-MEMS (micro- 
electronical-mechanical-systems) and bio-NEMS 
(nano-electronical-mechanical-systems) are attracting 
growing interest due to their huge potential in bio- and 
nanotechnology, e.g., for analyzing and manipulating 
tiny samples. Due to the large surface-to-volume ratios 
in nanoconfined systems, the flow in such devices is 
strongly influenced by the specific properties of the 
boundaries, i.e., by surface characteristics like the 
wetting behavior and/or slippage. 

If electric fields are involved, one particularly impor- 
tant mechanism is electroosmotic transport: In contact 
with a liquid, many materials commonly used in nan- 
otechnology (e.g., polydimethylsiloxane (PDMS)) ac- 
quire surface charges due to the ionization of surface 
groups [1]. Surfaces are thus often covered by a com- 
pensating counterion layer If one applies an ex- 
ternal electric field, the ions are driven in one direc- 
tion and drag the surrounding fluid along, thus creating 
the so-called electroosmotic flow (EOF). This electroki- 
netic effect has numerous consequences: It alters drasti- 
cally the migration dynamics of mesoscopic objects like 
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polyelectrolytes or colloids [3]. In microchannels, the 
EOF generated at the channel walls induces a total net 
flow, which is technologically attractive because it can 
be controlled and manipulated more easily on the sub- 
micrometer scale than pressure- or shear-driven flow. 

One important application of microchannels is to 
use them for separating different fragments of biolog- 
ical molecules like DNA by their length for sequenc- 
ing or further manipulation. In free solution, the elec- 
trophoretic mobility of high molecular weight polyelec- 
trolytes is length independent |3]. Electrophoretic sep- 
aration methods therefore often introduce the samples 
into micro- or nanostructured environments, e.g., dis- 
ordered gels (in sel electrophoresis), or structured mi- 
crochannels iSlillllSllililllillillil. The 

migration behavior of the molecules in such a setup re- 
sults from a complex interplay of electrostatics, hydro- 
dynamics, and confinement effects on the molecules, 
and modeling all of this in full detail is computation- 
ally challenging. Fortunately, it turns out that due to a 
combination of subtle screening effects - which are not 
yet fully understood - simple implicit-solvent Brown- 
ian Dynamics simulations that altogether neglect elec- 
trostatic and hydrodynamic interactions give results that 
are in qualitative and semiquantitative agreement with 
experiments ||3, 14, 15]. Nevertheless, such simplified 
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treatments miss important and interesting physics, and 
simulations with explicit solvent and explicit charges 
are clearly desirable. 

Here we present such explicit simulations of EOF 
and polyelectrolyte electrophoresis in nanochannels 
1 16, 17 1, using Dissipative Particle Dynamics (DPD) 
1 18l I19II . which is a popular mesoscopic simulation 
method. Particular emphasis is put on the role of the 
channel boundaries, i.e., on the effect of slippage and 
the electrical double layer at the surface. 

The remainder of the paper is organized as follows: 
After some theoretical considerations, we describe the 
simulation model and method, with special focus on our 
method to implement variable hydrodynamic bound- 
aries. Then we discuss the results, first for the sim- 
pler situation of a charged channel which only con- 
tains counterions, and finally for the full problem (elec- 
trolyte of varying ionic strength plus polyelectrolyte 
plus counterions). The results are compared with theory 
whereever possible. We conclude with a brief summary. 

2. Theoretical considerations: EOF in slit channels 

We consider a planar slit channel with identical walls 
at z = ±i/2, exposed to an external electric field in 
the X direction. The electrostatic potential O then takes 
the general form ^{x, y, z) - i//iz) + E^x + const, where 
we can set il/{0) - for simplicity. Comparing the 
Poisson equation for the electrostatic potential tfr with 
the Stokes equation and exploiting the symmetry of the 
channel, one finds IIT)] 



(1) 



where e,. is the dielectric constant, 77.5 the shear viscosity 
of the fluid, and v^^p = v^iO) is an integration constant. 

On the nanoscale, the appropriate hydrodynamic 
boundary condition at the channel walls is the partial- 
slip boundary condition 



6b d,v(x)\,„ = vAz)\z„ 



(2) 



where v^iz) denotes the component of the velocity in 
x-direction evaluated at the position zb of an effective 
"hydrodynamic boundary" position, and the second ef- 
fective parameter, the slip length Sb, characterizes the 
amount of slippage at the surface. Inserting Eq. ^ into 
Eq. ([T]i, we finally obtain the following simple expres- 
sion for the electroosmotic mobility. 



/^EOF = ^EoJE,, = ^p^p {l+K Sb), 



(3) 



where we have defined an inverse 'surface screening 



length' HZD 



and jJ^^^ is the well-known Smoluchowski result 
for the electroosmotic mobility at sticky walls fJ^^^ - 
—er 4'izB)hs- A similar result can be derived within 
the linearized Debye-Hiickel theory with k - kd (the 
Debye-Hiickel lengfli) Q. 

3. Simulation Method 

We study fluids that are driven through planar slit 
channels in the x-direction by external fields, applying 
periodic boundaries in the x and y dimension and re- 
pulsive (charged) walls in the z direction. The systems 
contained varying amounts of charged ions (anions and 
cations), plus possibly a polyelectrolyte chain. They 
were thermalized by the momentum conserving DPD 
thermostat 1. 1 8„ J,9], where the forces acting on parti- 
cles / are given by a sum = F'r + + Ffj) 
of standard conservative contributions Ff, dissipative 
forces FPj, and random forces Ff. with 



Fn = -JDPD (^(nj) (hj ■ Vij) fij 



(5) 
(6) 



K := +dz\pl\p\z=±ZB^ 



(4) 



Here aj(r) is an arbitrary weight function with finite 
range rc (chosen linear in our case, <jj{r) = 1 - r/r^ 
for r < rc), Jdpd is a friction coefficient, r,j = 
the unit vector in the direction of particle j, T the tem- 
perature, kB the Boltzmann factor, and (ij = ^j, are un- 
corrected Gaussian distributed random variables with 
zero mean and unit variance. In DPD simulations, the 
conservative forces are often taken to have a certain soft 
shape. Here we only use the DPD thermostat as de- 
scribed above. All simulations have been carried out 
with extensions of the freely available software package 
ESPResSo iiH. 

3.1. Simulation Model 

All particles, solvent, ions, and chain monomers, are 
modeled explicitly, and have the same mass m for sim- 
plicity. Ions and monomers repel each other with a soft 
repulsive Weeks-Chandler-Anderson (WCA) potential 
1 2211 of range o" and amplitude e. The same poten- 
tial acts between particles and the walls. In addition, 
chain monomers are connected by harmonic springs 
Uharmonic — jKhj - '"o)^ with the Spring coustaut k - 
25e/(T^ and ro = l.Oo". Neutral solvent particles have 
no conservative interactions except with the walls. The 
wall contains immobilized, negatively charged particles 
at random positions. Every second monomer on the 
polyelectrolyte carries a negative charge. All charges 
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Figure 1: Slip length rjg in units of 7,c vs. a for varying values of 
the parameter triplet (p, yopD, yt) (in units of tr^^ or 'Jmejcr, respec- 
tively). Black: series with p fixed: (3.75 , 2-10, 0.1-1). Red: series 
with jDPD fixed: (3.75-12.5, 2, 0.1-1). Blue: selected triplet values: 
(6.35,5,1),(5,5,1), (11.25,2,1.1), (11.25,2,1.2), (3.75,10,2.5). Dashed 
black line: Theory (Eq. 1101 ) The inset shows a blowup of the same 
data. After Ref [2311. 



are monovalent, and the system as a whole is elec- 
troneutral. Charged particles interact with each other 
via a Coulomb potential with the Bjerrum length Ab - 
lAnerksT - l.Ocr, and also with an external electric 
field - -l.Oe/ecr. Specifically, we show here results 
for systems with channel width 8cr and a surface charge 
density cta - -0.208ecr"^, which is found to correspond 
to the 'weak-coupling regime' il6ll . i.e., the regime 
where the Poisson-Boltzmann theory is valid. The to- 
tal counterion density was roughly p„umer = 0.06cr"^ and 
the salt density varied between p, =0.05625, 0.0375, 
0.03, 0.025, and 0.015cr"''. In molar units, this corre- 
sponds to 0.272, 0.181, 0.145, 0.121 and 0.072 mol/1, 
if we identify Ab ~ 0.7 nm, i.e., the Bjerrum length in 
water at room temperature [3J . 

3.2. Tunable Slip boundaries 

To realize arbitrary hydrodynamic boundary condi- 
tion at the walls, we introduce an additional coordinate- 
dependent viscous force that mimicks a wall/fluid fric- 
tion m 



F]- ^F? + Ff 



with a dissipative contribution 



(7) 



(8) 



coupling to the relative velocity (vi-v„aii) of the particle 
with respect to the wall, and a stochastic force 



(9) 



Here /? runs over/3 - x, y, z, Xifi is a Gaussian distributed 
random variable with mean zero and unit variance, and 




Figure 2: Flow profiles for counterion-induced flow at field strengths 
£v = 0.8 - l.OksT/e<T for varying slip lengths in the weak cou- 
pling regime. The straight lines represent the theoretical prediction 
ofEan.fTll. 



the weight function oJiiz) is chosen cjl(z) - I - zlzc 
up to a cut-off distance Zc- The prefactor ji sets the 
strength of the friction force and hence determines the 
value of the slip length. With this approach it is possible 
to tune the slip length 6b systematically from full-slip 
to no-slip and even (small) negative slip. Furthermore, 
one can use the Stokes equation to derive an analytical 
expression for the slip length 5^ as a function of the 
model parameters, giving ||23 1 



6_B 

z, 



— = -1 + 



(10) 



where F is the Gamma-Function, / the modified Bessel 
function of the first kind, and the dimensionless parame- 
ter a - zIjlP/Is depends on the density p and the shear 
viscosity tJs of the fluid. To test this expression, we have 
studied driven neutral fluids in our slit channels for the 
range of parameters -yi = (0.1 -5) ■\/rne/o-, p = (3.75- 
12.5) cr"^, and yoPD = (I-IO) ^/m€/c^. By perform- 
ing Plane Poiseuille and Plane Couette flow simulations, 
one can determine 77,, 6b, and the position of the hydro- 
dynamic boundary independently [23]. As expected, the 
hydrodynamic boundary is always close to the physi- 
cal boundary, whereas the slip length varies over a wide 
range. Fig. [1] shows that the results agree very nicely 
with the theoretical prediction, Eq. (fTOt lEsll . 

4. Results 

4.1. Counterion-induced electroosmotic flow 

We begin with discussing the "simpler" situation 
without polyelectrolyte and salt, where the fluid only 
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contains the counterions of the charges in the wall. In 
the regime of validity of the Poisson-Boltzmann equa- 
tion (the weak coupling limit) the potential distribution 
can be calculated analytically, giving 111] il/{z) oc 
log(cos^ (kcZ)), which results in the counterion distri- 
bution Pciz) = po/ cos^{KcZ) with the screening con- 
stant K-^ - e^po/lerksT . The parameter po (the coun- 
terion density in the middle of the channel) is set by the 
electroneutrality requirement, i.e., the integrated coun- 
terion density jdz pdz) must equal the surface charge 
density cr^. In our case, the corresponding calculation 
gives Po - 0.0174cr"-', in accordance with the numeri- 
cal value, po = (0.0176 + 0.0001)0--^ (l^. By virtue 
of Eq. ([TJ combined with one finally obtains the ex- 
plicit expression 

Vxiz) = . , „ EA log! — -) 

A-ttAeZtIs V ^ cos^(a-cZ) ' 

+2 /Cc(5b tan(AfcZB) j (11) 

for the flow profile, where we have expressed in terms 
of the Bjerrum length As- Fig. |2] compares our nu- 
merical results for varying surface characteristics (slip 
lengths) and field amplitudes with the theoretical pre- 
diction, Eq. (fTTT l. The theory describes the data very 
nicely, without any fit parameter. 

4.2. The full problem 

We are now ready to consider the full problem, i.e., 
a system containing polyelectrolyte, counterions, and 
varying amounts of salt ions, p, - 0.015 - 0.056cr"-'. 




Figure 3: Inset: Distribution of cations (salt cations and coun- 
terions, blue symbols) and anions p„ (salt anions, red symbols) for 
a system containing polyelectrolyte at the salt concentration p, = 
0.05625o""^. Main frame: Corresponding ionic difference profile 
Ap = pc - Pa The exponential fit (black line), gives the eft'ective in- 
verse screening length k = 2.305 ± 0.025(T"'. 
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Figure 4: Ratio l^^or /l^o eof plotted against Sbk for salt concentra- 
tions p, ranging in 0.015 0.056(T"^ and vaiious slip lengths. The 

blue line is the theoretical prediction of Eq. (|3] with slope 1 + SgK. 
Inset: Surface potential i^(zfi) as a function of salt concentration ps . 

Fig. |3] shows the ion distribution profiles for one salt 
concentration. The profiles of the ionic difference ex- 
hibit an exponential behavior f{z) = A(e^'^~ + e"^ + c), 
although the fitted screening parameter k - (2.305 + 
0.025)cr ' does not agree well with the Debye-Hiickel 
screening length, ko - 1.21cr"'. Similar observations 
were made at all other salt concentrations: For the 
surface charge ctj = 0.208cr"^, the linearized Debye- 
Hiickel theory is not valid. Nevertheless, a well-defined 
surface screening length k can be extracted from the data 
by a simple exponential fit. 

Next we discuss the electroosmotic flow in these 
channels. Fig. |4] compiles our numerical results for 
the EOF mobility for all salt concentrations and slip 
lengths. They are in very good agreement with the the- 
oretical prediction of Eq. ([Sj, where has been de- 
termined independently by a linear regression for each 
salt concentration. It is worth noting that the presence 
of the polyelectrolyte does not perturb the amplitude of 
the electroosmotic flow. 

Finally, we consider the effective migration of the 
polyelectrolyte in the electric field. It results from a 
combination of two eff'ects: the 'bare' electrophore- 
sis relative to the surrounding fluid, and the convective 
transport by the electroosmotic flow. In many situations 
of interest, one can arg ue that these two contributions 
simply add up ifTil [l7n . The total mobility p, is then 
given by /i, = ^^^^ + fi^, where p^ is the electrophoretic 
mobility of the polyelectrolyte in a hypothetical fluid 
at rest, and by virtue of Eq. (jS), it can be expressed in 
terms of the electroosmotic mobility p^^p as 

J^^l + -—^ , (12) 

/^OF PeOf(1 
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Figure 5: Ratio //f/|/jj;Qf.| plotted against SgK for all salt concentra- 
tions. The black line is the theoretical prediction of Eq. | |121 . with 
one single fit parameter /"e/p^j,^ = -3.778 ± 0.128. Negative values 
of /I, /ll^^Qp I indicate absolute negative total mobilities of the poly- 
electrolyte. Inset: Total displacement of the polyelectrolytes center 
of mass for different boundary conditions at the salt concentration 
p, = 0.05625o""'. The lines coiTespond from top to bottom to the 
slip lengths Sb » (0.00, 1.292, 1.765, 2.626, 5.664, 14.98)d-. 



where the ratio jj^/fi'^^^ is expected to depend only 
weakly on the ionic strength of the electrolyte and the 
slip length of the surface. The main effect of slippage is 
incorporated in the factor (1 -i- a- 5b) ' OlTll . 

Our numerical results for the total mobility of the 
polyelectrolyte for varying boundary conditions are pre- 
sented in Fig. |5] They are in excellent agreement with 
the theoretical prediction, Eq. (fTSl i. with one single fit- 
ted ratio /"f//""of = -3.778 + 0.128. For stick bound- 
aries {6b ~ 0) one obtains ordinary behaviour where the 
polyelectrolyte follows the electric force acting on the 
monomers. In the presence of slip, the absolute mo- 
bility may become negative if the electroosmotic flow 
exceeds a critical value. This is because the immobile 
wall charges and the charges on the polyelectrolyte have 
the same sign, hence the directions of the EOF and the 
bare electrophoresis are opposite. If the wall charges 
and the polyelectrolyte charges are opposite, slippage 
effects should enhance the total mobility of the poly- 
electrolyte. 



5. Conclusions and Outlook 

We have presented mesoscopic DPD simulations of 
EOF and polyelectrolyte electrophoresis in narrow mi- 
crochannels, taking full account of hydrodynamic and 
electrostatic interactions. Slippage effects massively in- 
fluence the electroosmotic flow and therefore the to- 
tal mobility of the polyelectrolyte. Under certain con- 
ditions, even a negative mobility can be achieved, in 



agreement with recent experiments OlSll . All our nu- 
merical results are in good agreement with analytical 
expressions, which were derived based on the Stokes 
equation. 

In sum, our mesoscopic simulations indicate that the 
migration of polyelectrolytes in nanochannels results 
from the interplay of electroosmotic, electrophoretic, 
electrostatic and slippage effects. To describe the mobil- 
ity adequately, all of these factors need to be accounted 
for. From a technological point, the characteristics of 
the channel walls could be used to significantly enhance 
flow amplitudes, which offers the possibility to reduce 
the time needed for polymer migration or separation 
techniques. This could be an important aspect for fu- 
ture applications in microchannels or micropumps. 
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